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ABSTRACT 

The  probability  of  m  correlated  random  variables  drawn  from  a  multivariate 
normal  distribution  being  non-negative  is: 

[CO  fOO  fOO 

Jo  Jo  Jo  . *n)3*i,d*2 . dxn. 

Exact  results  for  this  probability  integral  are  unavailable  for  m  >  3.  Ap¬ 
proximations  for  higher  dimensional  problems  have  generally  yielded  poor 
results  except  for  special  cases,  such  as  compound  symmetry,  which  is  of 
limited  value  in  practice.  The  purpose  of  this  paper  is  to  present  a  gen¬ 
eral  approximation  of  this  probability  integral.  The  algorithm  developed 
here  is  computationally  tractable  for  m  =  50  and  accurate  for  very  general 
correlational  structures.  The  performance  of  this  algorithm  is  compared  to 
results  based  on  Clark’s  (1961)  original  approximation,  Gaussian  quadrature 
formulae,  and  Monte  Carlo  simulation  methods.  Application  of  this  approx¬ 
imation  to  problems  of  conditional  dependence  in  IRT  estimation  problems 
is  discussed. 


1  INTRODUCTION 


With  the  renewed  interest  in  marginal  maximum  likelihood  methods  (MML) 
in  item  response  theory  (Bock  and  Aitken,  1981),  the  necessity  for  approx¬ 
imating  orthant  probabilities  of  the  multivariate  normal  distribution  has 
arisen.  In  the  MML  estimation  of  item  parameters,  for  example,  all  *  <  2n 
orthant  probabilities  must  be  evaluated  in  order  to  obtain  the  likelihood 
of  the  parameters  of  a  given  IRT  model.  Similarly,  the  same  multivariate 
normal  orthant  probabilities  are  required  to  obtain  Bayesian  estimates  of 
ability  assuming  dichotomous  scoring. 

Because  IRT  models  assume  conditional  independence  among  the  items, 
that  is,  all  association  among  the  n  items  is  completely  explained  by  their 
joint  association  with  the  latent  variables,  m-dimensional  Gaussian  quadra¬ 
ture  formulae  have  been  used  to  provide  the  necessary  orthant  probability 
estimates.  In  this  case,  the  simple  “product  formulae”  for  numerical  inte¬ 
gration  can  be  applied,  because  the  residuals  of  the  model  have  the  simple 
uncorrelated  multivariate  normal  distribution  e  ~  (0 ,<r?J),  where  o*  is  the 
so-called  item  “uniqueness” ,  that  is,  1  -  a*,  and  o tj  is  the  item  factor  loading. 
However,  these  numerical  integrations  become  extremely  expensive  as  the 
dimensionality  of  the  problem  increases.  To  the  extent  that  the  items  have  a 
simple  factor  structure,  these  approximations  are  accurate  to  any  practical 
degree,  however  their  behavior  with  more  general  correlational  structures  is 
predictably  poor  and  the  extent  of  bias  introduced  into  parameter  estimates 
is  unknown. 


In  the  general  case,  the  probability  that  m  correlated  random  variables 
drawn  from  a  multivariate  normal  distribution  are  non-negative  is: 


r  r...r 

Jo  Jo  Jo 


(*i,  *J,  •  •  •  ,*m)  dzu dx2, . . . ,  dz„ 


Exact  results  for  the  bivariate  normal  distribution  have  been  obtained  from 
the  work  of  Sheppard  (see  Kendall  and  Stuart,  1973).  David  (1953)  has 
shown  that  trivariate  normal  probabilities  are  a  direct  extension  of  Shep¬ 
pard’s  earlier  result.  However,  exact  results  for  m-variate  normal  probabil¬ 
ities  for  m  >  3  are  unknown. 


Kendall  (1941),  McFadden  (1960)  and  Moran  (1948)  have  developed  in¬ 
finite  expansions  of  the  quadrivariate  normal  integral,  but  these  results  are 
computationally  cumbersome.  More  tractable  representations  have  been 
proposed  by  Abrahamson  (1964),  Dutt  (1973),  Dutt  and  Lin  (1975)  and 
Childs  (1967).  Even  in  the  quadrivariate  case,  however,  these  approxima¬ 
tions  are  intractable  for  general  correlational  structures  and  have,  therefore, 


only  been  applied  under  restrictions  of  compound  symmetry  or  band  matri¬ 
ces  in  which  p  =  .5  just  off  the  main  diagonal. 

In  a  previous  paper,  we  have  developed  an  IRT-type  model  for  estimating 
trend  in  correlated  proportions  (Gibbons  and  Bock,  1987).  In  this  case  the 
binary  items  are  not  components  of  the  same  test,  but  rather,  dichotomous 
classifications  orde?  _d  in  time.  Since  the  association  between  temporally 
proximal  responses  will  typically  be  greater  than  the  association  between 
temporally  distal  responses,  the  assumption  of  conditional  independence  is 
untenable  and  residual  autocorrelation  must  be  assumed.  At  the  sugges¬ 
tion  of  James  Heckman  (personal  communication)  we  adapted  the  so-called 
Clark  Algorithm  (Clark,  1961)  to  the  problem  of  evaluating  the  likelihood  of 
this  model.  Our  modified  Clark  Algorithm  combines  both  Gaussian  quadra¬ 
ture  formulae  with  Clark’s  original  approximation  which  is  extremely  fast. 
The  purpose  of  this  paper  is  to  present  details  of  this  modified  approxima¬ 
tion,  illustrate  its  computation,  and  determine  its  accuracy  in  a  few  relevant 
examples. 

2  THE  CLARK  ALGORITHM 

Designating  positive  directions  1  and  negative  directions  0,  we  may  repre¬ 
sent  the  probability  of  the  positive  orthant  of  an  m-variate  distribution  by 
P(l,  1, ....  1),  that  of  the  negative  orthant  by  P(0, 0, . . . ,  0),  and  that  of  any 
one  of  the  other  2m  -  2  orthants  by  inserting  the  appropriate  pattern  of  l’s 
and  0’s.  The  Clark  algorithm  provides  a  computing  approximate  for  any 
orthant  of  a  multivariate  normal  distribution  with  arbitrary  vector  mean 
and  covariance  matrix.  Clark  (1961)  derives  the  following  formulas. 

Let  any  three  successive  components  from  an  m-variate  vector,  |/t ,  be 
distributed: 

Vi  /[  * 

V.+  l  ~  N  I  A*i+ 1  i  <7\&i+lPi,t+l  &t+l 

Vi+2  V  + 2  .  a*a*+2Pi.i+2  ffi+lffi+2Pi+l,i+2  <*i+ J 

Let  y*  =  max(y,)  =  y»,  and  compute  the  probability  that  yj+i  >  y<  as 
follows: 


I 

1 


set 

*i+ 1 

=  (Mi  ~  Mi+ i)/ft+i, 

where 

&  i 

=  fff  +  <r?+1  -  2<r<<r<+1p,i<+1 

Then 

P(t*«  >  if) 

=  ^(y<+i  -  y  >  0) 

=  *(-«♦!) 

the  value  of  the  univariate  normal  distribution  function  at  the  standard 
deviate  -*j+ 1. 

Now  let  ii+i  =  max(]ft,y,>i)  and  assume  (as  an  approximation)  that 
(v»+J»  Vi+i)  »  bivariate  normal  with  means, 

m(w+i)  =  £(Vi+j)  =  Mi+i 

M(ii+ 1)  =  £($i+i)  =  W*(a+i)  +  P»+i*(-*<+i) 

variances 


*-2(y.+j)  =  £(vhi )  ~  tHtfi+l)  =  *i+2> 
*2(y.+i)  =  £(y?+1)  -  f2(y.+i), 

where 


£  (y.>i)  =  (m!  +  +  (m*+i  +  *2+i)$(-*<+t)  +  (Mi  +  Mi+l)Si+i<f>(zi+i), 

and  correlation 


»(y.+i) 


p(y<+i>y.+2)  = 

Then, 

P(y,+a  =  max(y,,y<+i,y<+j))  =  P((y,+j  -  y,+i  >  0)  n  (y,+2  -  Vi  >  0)) 
is  approximated  by 


•P(y<+2  >  y»+i)  =  P(y.+:  -  y«+i  >  o) 


Mi+2  -  M(i/i+ 1) 


V**>»  +  **(fc+l)  -  2<r,+j<r(y,+1)p(y<+1,y,+2) 


$ 


Assuming  as  a  working  approximation  that  &+i  is  normally  distributed 
with  the  above  mean  and  variance,  we  may  therefore  proceed,  recursively 
from  *  =  1  to  *  =  m  —  1,  where  ym+i  is  *n  independent  dummy  variate  with 
mean  zero  and  variance  zero  (i.e.  ym+i  =  0). 

Then 


^(Vm+i  =  max(yi,yj, .  •  -  <  Vm+x)) 

=  />((Vm+l -yi  >  0)n(ym+1  -  y,  >  0)n...n(ym+1  -  ym  >  o)) 
=  P({-vi  >0)n(-yj  >0)n...n(-ym  >  o)) 

approximates  the  probability  of  the  negative  orthant  of  the  specified  mul¬ 
tivariate  normal  distribution.  The  probability  of  any  other  orthant  can  be 
obtained  by  reversing  the  signs  of  the  variates  corresponding  to  l’s  in  the 
orthant  pattern. 

3  COMPUTATIONAL  EXAMPLE 

For  computational  purposes,  it  is  convenient  to  set  $(-z)  =  1  -  $(z)  and 
rewrite  Clark’s  equations  as: 


M(SN+l)  =  Mi+x  +  (Mt-/i<+l)$(*+i)  +  ?,+i<£(*i+i), 

£(y*+i)  =  M?+i  +  ehi  +  (m?  +  e*  ~  M?+i  -  <r»>i)$(z,+i) 

+0*  +Mi+lki+l^(*.+i). 


a  (y.+i)  —  ^(y*+i)  -  f*(y.+i), 


^J(y<+x,yi+j) 


^J(y.+i.y.+j)  +  kJ(fc,y<+2) 

-<7J(y.+i,yt+2)]*(si+i) 


Given  this  transformation,  the  ith  step  now  only  requires  a  single  evaluation 
of  $(.)  and  <£(.).  For  example,  suppose  that 


M  = 


1 


r 


3 

2 

2 


and  E  = 


3 

1  2 
1  0  2 


Given  these  values  for  fi  and  E,  the  probability  of  the  negative  orthant  (-  - 
•)  can  be  obtained  as  follows. 

First  find  P(y s  >  yj)  as: 


ft 

z2  = 

£(&)  = 
^ (vl)  = 

<r2(h) 

Therefore; 

P{ ys  >  yj) 


V^l  +  °\  ~  2<t13 
^3  +  2  -2(1)  =  1.732 

(mi  -  M2)/ft 

(3  -  2)/1.732  =  .577 

M2  +  (mi  -  M2)$(*2)  +  f2<£(*2) 

2  +  (3  -  2)$(.577)  +  1.732fli(.577)  =  3.30 


Mj  +  o\  +  (m?  +  -  M2  ~  *1  )$(*2) 

+  (Ml  +M2)?2^(*2) 

6  +  6$(.577)  +  5(1.732)*(.577)  =  13.23 


£(v  i)  -  **(M 

13.23  -  3.302  =  2.34 


*23  +  (*I3  -  *fs)*(*2) 

$(.577)  =  .72 


=  $ 

=  $ 


2-3.3  \ 

\/2  +  2.34  -  2(.72) ) 


.222 


To  determine  P(y j  >  0)  we  set  y<  =  0  and  <r2  =  0,  and  compute: 


*s  =  (m 2  -  Ms) /fs 

=  (3.3  -  2)/1.703  =  .763 

^  (to)  =  MS  +  (#2  -  Ms)$(*s)  +  ft<K*s) 

=  2  +  (3.3  -  2)$(.763)  +  1.703<£(.763)  =  3.518 

t (Vs)  =  4  +  4  ’  (Mj  +  -  4  ~  ffs)$(*s) 

+(£2  +  Msjft^s) 

=  6  +  7.23$(.763)  +  5.3(1.703)^(.763)  =  14.307 

**(»)  =  ^(y|)  -  ^2(ys) 

=  14.307  -  3.518*  =  1.931 

<TJ(ys,y4)  =  <Tu  +  (*14  -  ^L)*(*s) 

=  0 

Therefore; 


^(y4  >  y3)  = 


M<  ~  M3 

/<r*+H  ~  2fflt 


-  J  °~3  518 _ w*>‘ 

"  \y0+ 1.931-2(0)]  ~  1 


Hence,  P(0  >  yj)  =  $(-2.53)  =  .006,  which  is  the  probability  of  orthant  (- 
-*)• 


4  THE  MODIFIED  CLARK  ALGORITHM  #1 

In  our  previous  paper  (Gibbons  and  Bock,  1987),  we  noted  that  the  accu¬ 
racy  of  the  Clark  approximation  diminishes  with  increasing  magnitude  of 


the  correlations.  If  we  apply  the  Clark  approximation  directly  to  estimates 
of  inter-item  correlations,  it  will  generally  yield  biased  results  due  to  the 
size  of  correlations.  This  is  true  regardless  of  whether  the  correlation  ma¬ 
trix  exhibits  the  property  of  conditional  independence.  Alternatively,  if  we 
examine  the  residual  inter-item  correlation  matrix  at  fixed  points  on  the 
ability  scale,  we  will  observe  the  identity  matrix  for  conditionally  indepen¬ 
dent  solutions  or  small  residual  correlations  for  those  item  pairs  that  are 
conditionally  dependent.  In  light  of  this,  we  evaluate  the  response  function 
at  several  fixed  points  on  the  ability  scale  using  Gauss-Hermite  quadrature, 
and  correct  these  estimates  using  the  Clark  algorithm.  These  corrections  de¬ 
pend  only  on  the  residual  inter-item  correlations,  which  in  practice  should 
be  quite  small.  The  first  modified  Clark  algorithm  proceeds  as  follows. 

Step  1  Obtain  a  factor  solution  of  dimension  m,  using  full  information  fac¬ 
tor  analysis  for  binary  data  (Bock  and  Aitken,  1981,  and  Bock,  Gibbons 
and  Muraki,  1986)  if  the  correlations  are  unknown  or  using  principal  factor 
analysis  if  the  item-correlations  are  known,  as  in  the  following  simulations. 

Step  2  Using  the  estimated  factor  loading  matrix  Anxk  compute  the  esti¬ 
mated  residual  correlation  matrix  /%xn  as  At,Ay  for  »  ^  j  else 

P'ij  =  1>  where  t  is  the  dimensionality  of  the  space  we  are  conditioning  on 
plus  1.  For  example,  in  the  unidimensional  case,  t  =  2.  This  correlation 
matrix  represents  the  degree  of  residual  conditional  dependence. 

Step  S  Given  the  previous  values  of  item  thresholds  7,  and  item  factor  load¬ 
ings  A ,y  for  the  t  -  1  prinical  factors,  compute  the  invariant  item  parameters 
dj  (slope)  and  6y  (intercept).  For  example,  in  the  unidimensional  case, 

a,  =  VyTTx? 
and 

b)  =  -V\A-a> 

Step  5  At  each  point  on  the  ability  dimension  (i.e.  at  each  quadrature  node 
Xk  )  compute  the  value  of  the  response  function  for  each  item  as: 

z]k  =  Cj  +  GjXk 

where  Cj  =  -ajbj  and  Xk  are  the  nodes  of  the  Gauss-Hermite  polynomial 
(see  Stroud  and  Sechrest,  1966). 


Step  6  At  each  quadrature  point,  substitute  the  values  of  zy*  for  the  mean  | 

vector  (a  and  R'  for  the  covariance  matrix  E  and  compute  the  Clark  approxi-  j 

mated  probability  Cj(X*).  Accumulating  these  probabilities  over  all  quadra-  i 

ture  nodes  for  a  given  response  pattern  (x*)  yields  the  desired  marginal 
probability  estimate 

j 

I 

Ma)  =  r  <hm*m  ' 

J -00  • 

=  Ec,(**m(**)  j 

*- 1  i 

where  A(Xk)  is  the  corresponding  weight  at  quadrature  node  A*.  J 

We  note  that  in  practice,  the  effect  of  assuming  normality  of  the  max-  ' 

imum  of  two  jointly  normal  variables,  produces  probability  overestimates 
in  the  tails  of  the  distribution.  As  such,  we  apply  an  empirically  based 
correction  factor  to  these  probability  estimates  which  involves  raising  the 
Clark  adjusted  probability  estimate  to  the  power  1.3.  This  correction  factor 
appears  to  provide  the  neccessary  adjustment  across  the  entire  quadrature 
space. 

5  THE  MODIFIED  CLARK  ALGORITHM  #2 

Algorithm  #1  is  computationally  expensive,  because  the  Clark  approxima¬ 
tion  must  be  evaluated  for  each  response  pattern  at  every  point  in  the 
quadrature  space.  In  a  problem  with  1000  unique  response  patterns  and 
two  principal  factors,  each  with  10  quadrature  points,  the  Clark  approxima¬ 
tion  must  be  invoked  100,000  times  per  iteration.  An  alternate  approach  is 
to  simply  use  the  Clark  approximated  probability  as  a  correction  term,  ap¬ 
plied  directly  to  the  usual  probability  estimate  obtained  from  the  quadrature 
solution;  that  is: 
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Returning  to  the  previous  example  of  1000  unique  response  patterns  and 
10  quadrature  points  in  each  of  two  dimensions  (i.e.  100  points  in  total),  the 
Clark  algorithm  need  only  be  invoked  1000  times  per  iteration;  that  is,  once 
for  each  response  pattern  in  contrast  to  100,000  times  for  algorithm  #1. 

The  probabilities  obtained  using  algorithm  #2  differ  from  thoee  obtained 
from  algorithm  #1  in  that  they  will  not  sum  to  unity  even  if  all  response 
patterns  are  realized  in  the  sample.  When  the  number  of  items  is  small,  say 
ten  or  less,  these  probabilities  can  be  normalized  to  yield  the  appropriate 
metric;  however,  in  larger  problems  (i.e.  with  eleven  or  more  items)  normal¬ 
ization  is  not  possible,  because  all  patterns  generally  are  not  realized  in  even 
large  samples.  We  note,  however,  that  maximum  likelihood  estimation  does 
not  require  normalized  probabilities  as  long  as  their  relative  magnitudes  are 
invariant  to  transformation  of  scale.  This  condition  does  hold  for  algorithm 
#2. 

6  ILLUSTRATION 

To  examine  the  accuracy  of  the  modified  Clark  Algorithm,  we  designed  the 
following  limited  simulation  study.  First,  we  simulated  one  million  five- 
variate  normal  deviates  for  each  of  the  following  conditions: 

1)  compound  symmetric  matrices  with  p  =  0.2  through  p  =  0.8.  For  exam¬ 
ple. 


'  1.0 

0.5 

1.0 

0.5 

0.5 

1.0 

0.5 

0.5 

0.5 

1.0 

0.5 

0.5 

0.5 

0.5 

t)  Autocorrelated  matrices  with  p  =  0.2  through  p  =  0.8.  For  example, 
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S)  Conditionally  dependent  correlation  matrices  with  principal  factor  load¬ 
ings  of  Au  =  An  =  Ais  =  Am  =  Am  ranging  from  0.5  to  0.7  and  two  method 
related  factors,  the  first  with  An  =  An  =  A«  ranging  from  0.2  through  0.7 
and  A]4  =  Am  =  0.0,  and  the  second  with  Aji  =  Ajj  =  Ajj  =  0.0  and 
Aji  =  A35  ranging  from  0.2  through  0.7. 

For  example,  the  correlation  matrix  corresponding  to  factor  pattern  ma¬ 
trix: 


0.7  0.3  0.0 
0.7  0.3  0.0 
0.7  0.3  0.0 
0.7  0.0  0.3 
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0.49  0.49  0.49  0.58 


In  each  of  the  above  simulated  conditions,  the  mean  vectors  were  zero, 
therefore,  the  binary  response  patterns  were  obtained  by  dichotomizing  the 
simulated  normal  deviates  at  zero.  The  dichotomized  response  patterns 
were  then  sorted  and  unique  patterns  and  their  respective  frequencies  were 
accumulated.  When  divided  by  one  million,  these  frequencies  yield  the  so 
called  Monte  Carlo  probability  estimates  which  should  be  exact  to  at  least 
four  decimal  places;  that  is,  the  standard  error  of  the  simulated  probabilities 
is: 


^p,(i -Pi)y 


which  has  a  maximum  value  of: 


1000000 


= .0005 


Accuracy  of  the  estimated  probabilities  was  determined  by  computing 
the  average  absolute  deviation  for  each  method  (i.e.  the  average  absolute 
difference  between  Monte  Carlo  and  Clark  estimates  over  all  response  pat¬ 
terns). 

In  an  effort  to  examine  the  properties  of  these  approximations  in  larger 
sets  of  items,  that  are  more  typical  in  practice,  we  produced  two  10-item 
simulations.  The  first,  consisted  of  ten  million  multivariate  normal  deviates 
for  the  compound  symmetric  case  (p  =  0.5),  and  the  second  for  the  condi¬ 
tionally  dependent  example  with  principal  factor  loadings  of  Ai,-  =  0.6  and 
method  related  factor  loadings  of  Aji  . . .  Aj.io  —  0.4  and  Aj.n . . .  Aj  jo  =  0.4. 

7  RESULTS 

Results  of  the  simulations  are  displayed  in  Tables  1-5.  Inspection  of  Table 
1,  which  displays  results  for  the  five-item  compound  symmetric  case,  reveals 
that  on  average,  both  the  quadrature  solution  and  the  modified  Clark  al¬ 
gorithm  #2,  recover  the  Monte  Carlo  estimates  to  four  decimal  places  (i.e. 
average  difference  p  =  0.5,  .0002)  whereas  the  modified  Clark  algorithm  #1 
is  slightly  less  accurate  (i.e.  average  difference  p  =  0.5,  .0009).  In  contrast, 
the  original  Clark  algorithm  is  inferior  to  the  other  methods  and  this  in¬ 
feriority  increases  with  larger  correlations  (p  =  0.8,  average  difference  = 
.0163).  These  results  are  exactly  as  expected  for  a  conditionally  indepen¬ 
dent  correlational  structure.  In  contrast,  inspection  of  Table  2  reveals  that 
autocorrelation  produces  inferior  estimates  for  the  quadrature  solution  (eg. 
p  =  0.5,  average  difference  =  .0070)  whereas  the  modified  Clark  algorithms 
produced  more  consistent  estimates  across  the  entire  range  (eg.  p  =  0.5, 
average  difference  for  algorithm  #1  =  .0041  and  for  algorithm  #2  =  .0036). 
Algorithm  #2  was  in  general,  slightly  better  than  algorithm  #1.  The  orig¬ 
inal  Clark  algorithm  produced  reasonable  estimates  through  p  =  0.5,  but 
deteriorated  quickly  for  values  of  p  >  0.5.  Since  the  elements  of  the  correla¬ 
tion  matrix  become  more  homogeneous  for  extreme  values  of  p  (eg.  p  —  0.8) 
it  is  not  surprising  that  the  performance  of  the  quadrature  solution  stabilized 
for  values  of  p  >  0.6. 

In  terms  of  five-item  conditionally  dependent  correlational  structures 
(see  Table  3),  the  modified  Clark  algorithms  performed  similarly  whereas 
the  original  Clark  algorithm  and  the  quadrature  solution  were  consistently 


t 


inferior,  for  even  moderate  dependence  (Aj;  or  As;  >  0.3).  The  modified 
Clark  algorithm  #2  produced  average  differences  that  were  slightly  smaller 
than  algorithm  #1  and  as  little  as  one-sixth  the  size  of  the  standard  quadra¬ 
ture  solution.  Overall,  performance  was  better  for  solutions  in  which  the 
principal  factor  loadings  were  smaller  (eg.  An  =  0.5,  An  =  0.5,  A31  =  0.0, 
average  difference  for  algorithm  #1  =  .0034,  average  difference  for  algorithm 
#2  =  .0031,  average  difference  for  original  Clark  =  .0041  and  average  differ¬ 
ence  for  quadrature  =  .0084  versus  An  =  0.7,  An  =  0.5,  A31  =  0.0,  algorithm 
#1  =  .0068,  algorithm  #2  =  .0048,  original  Clark  =  .0091,  and  quadrature 
=.0099). 

Results  for  the  larger  item-sets  are  presented  in  Table  4  for  the  compound 
symmetric  case  and  Table  5  for  the  conditional  dependent  case.  In  general, 
results  for  the  larger  item-sets  parallel  those  of  the  smaller  item-sets.  For 
the  compound  symmetric  case  (p  =  0.5)  the  average  differences  were  .000065 
for  algorithm  #1,  .000007  for  algorithm  #2,  .000007  for  the  quadrature 
solution.  For  the  conditionally  dependent  case  (Ay  =  0.6,  Ay  =  0.4  or  0.0 
and  A3;  =  0.4  or  0.0),  the  average  differences  were  .000218  for  algorithm 
#1,  .000242  for  algorithm  #2,  .000329  for  the  quadrature  solution. 

8  DISCUSSION 

The  results  of  this  study  clearly  demonstrate  that  when  the  assumption  of 
conditional  independence  is  violated,  bias  in  the  standard  IRT  probability 
estimates  are  produced.  Both  modified  Clark  algorithms  presented  here 
minimize  this  bias;  however,  algorithm  #2  produces  slightly  more  accurate 
results  than  algorithm  #1  for  small  numbers  of  items,  at  a  remarkable  sav¬ 
ings  in  computation.  Algorithm  #2  decreases  the  required  computation  by 
a  factor  of  qm ,  where  m  is  the  number  of  underlying  dimensions  and  q  is 
the  number  of  quadrature  points  in  each  dimension.  When  compared  to 
the  standard  quadrature  solution,  algorithm  #2  requires  an  additional  s 
evaluations  of  the  Clark  algorithm,  where  $  is  the  number  of  uniquely  ob¬ 
served  response  patterns.  Conversely,  when  method  related  factors  do  exist, 
multiple-factor  solutions  will  greatly  increase  the  computational  complex¬ 
ity  of  the  standard  IRT  approach,  whereas  no  increase  in  computation  is 
required  for  algorithm  #2. 

It  is  important  to  point  out  that  the  use  of  these  approximations  should 
not  replace  multiple  factor  solutions  where  the  additional  factors  contribute 
to  our  understanding  of  individual  differences.  Indeed,  the  methods  de¬ 
scribed  here  treat  these  potentially  meaningful  effects  as  errors  of  measure- 


ment.  By  allowing  for  more  general  typee  of  measurement  variability,  these 
methods  can  mask  important  characteristics  of  item-person  interactions. 
There  are,  however,  situations  where  small  method  related  effects  confound 
our  ability  to  accurately  characterize  the  dominant  dimension  of  interest, 
because  the  inter-item  association  is  not  strictly  a  function  of  a  single  un¬ 
derlying  trait.  In  these  cases,  the  ability  to  segment  these  more  complex 
measurement  errors  from  our  estimates  of  the  central  trait  or  aptitude  of 
interest  is  a  highly  desirable  goal  of  measurement. 

Our  future  research  in  the  application  of  the  Clark  algorithm  for  likeli¬ 
hood  evaluation  of  IRT  models,  will  focus  on  the  estimation  of  ability.  When 
method  related  factors  violate  the  assumption  of  conditional  independence, 
current  methods  for  estimating  ability  (Bock  and  Aitken,  1981)  will  pro¬ 
duce  incorrect  results.  Some  preliminary  work  in  this  area  suggests  that  the 
variance  in  Bayes  “expected  a  posterior?  estimates  (EAP)  for  fixed  levels 
of  ability,  increases  by  a  factor  of  2  to  3  in  the  presence  of  even  moderate 
dependence.  Substituting  Clark  estimates  for  the  conditional  probabilities 
that  are  usually  employed  in  calculating  EAP  estimates  should  account  for 
this  “extra”  normal  variability,  and  therefore,  provide  more  accurate  ability 
estimates. 

In  addition,  we  will  also  further  explore  the  difference  between  our  two 
modified  Clark  algorithms,  by  focusing  on  their  behavior  at  individual  points 
in  the  quadrature  space.  In  this  way,  we  hope  to  even  further  improve  the 
performance  of  algorithm  #1  relative  to  algorithm  #2,  and  perhaps,  develop 
a  third  procedure  which  is  an  improvement  over  both. 
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